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Recent studies on epistatic networks of model organisms have unveiled a certain type of modular 
property called monochromaticity in which the networks are clusterable into functional modules 
that interact with each other through the same type of epistasis. Here we propose and study 
three epistatic network models that are inspired by the Duplication-Divergence mechanism to gain 
insight into the evolutionary basis of monochromaticity and to test if it can be explained as the 
outcome of a neutral evolutionary hypothesis. We show that the epistatic networks formed by 
these stochastic evolutionary models have monochromaticity conflict distributions that are centered 
close to zero and are statistically significantly different from their randomized counterparts. In 
particular, the last model we propose yields a strictly monochromatic solution. Our results agree 
with the monochromaticity findings in real organisms and point toward the possible role of a neutral 
mechanism in the evolution of this phenomenon. 
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I. INTRODUCTION 

Highly interacting molecular networks (Tj-Q ^ire at the 
forefront of systems biology and are essential to un- 
derstand biological systems beyond the single molecule 
framework. In these networks usually a modular topol- 
ogy HI is predictive of the presence of functional mod- 
ules and therefore topological considerations [ol, [lo| have 
relevance in detecting the function and dynamics of these 
networks. The origin of the functional modularity is not 
thoroughly characterized: at first it was considered to be 
due to evolutionary pressure [7] but more recently obser- 
vations have shown that modularity might be the natural 
outcome of neutral evolution jTlj . 

Recently the attention has been directed to the study 
of complex genetic interactions [l2|, [l3| . In fact the rela- 
tionship between phenotypes and genotypes is a complex 
one - phenotypes are most commonly determined by si- 
multaneous interactions between multiple genes rather 
than by a single gene [l^. These genetic interactions 
are furthermore complicated by the fact that the effects 
of individual genes are modified by other genes, result- 
ing in the phenomenon called epistasis. The epistatic 
network between mutations or genetic variations in the 
genome is found to play a crucial role in determining the 
fitness of different organisms [l2|, |l3|. This network is 
interacting with the other biological networks of the cell 
such as the protein-protein interaction network [l5[ [l6| 
the metabolic network [it) or the transcription network 
[sl, UH and is crucial in determining the phenotype of 
an organism. The decree distribution of these networks 
is clearly fat-tailed [iSl thus showing a similar topology 
with respect to other biological networks. Moreover the 
genes involved in epistatic interactions can have a hierar- 
chical relationship; they can mask or modify each other's 
effects and even combine to create new phenotypes. 

A mathematical definition of epistasis that was intro- 
duced by population geneticists in order to deal with 
complex traits such as human diseases in a quantitative 
and statistical way is known as the multiplicative model. 



Widely used in genetic interaction studies, it makes use 
of the quantitative measure e which is the difference be- 
tween the fitness Wab of a double mutant and the prod- 
uct of fitnesses of two single mutants Wa-, Wb-, i.e. 

e = Wab - WaWb 

The sign of this deviation from multiplicative behavior 
helps distinguish between two different classes of interac- 
tions: a positive (antagonistic) epistasis signifies a buffer- 
ing type of interaction in which the effect of the double 
mutation is less severe than single mutations combined; 
a negative (synergistic) epistasis means an aggravating 
type of interaction where the double mutation results in 
a more severe effect on the fitness. 

Biological systems are known to often display a func- 
tionally modular architecture on various levels [l9|. 
Much of the robustness and adaptability of biological 
networks has been attributed to this modular structure. 
Modularity in biological networks has been studied from 
the perspectives of topology [205 , information theory [2l| 
and more recently, epistasis. 

A recent study [12| has investigated epistasis through 
a genetic interaction network constructed by the calcu- 
lation of epistasis values e from the fitnesses of all of the 
single and double knockouts of 890 metabolic genes of 
S. Cerevisiae (budding yeast). In both supervised orga- 
nization with known gene annotations and unsupervised 
organization with unknown gene functions, the study re- 
vealed a structural clusterability into modules that inter- 
act with each other monochromatically, that is, via the 
same type of epistasis. Monochromaticity in the inter- 
actions between functional modules was later verified as 
part of a genome-wide genetic interaction study [l3| of 
S. Cerevisiae. In this larger scale research, the genetic in- 
teraction map of the budding yeast was constructed from 
5.4 million gene pairs. Modules belonging to different 
pathways and complexes were shown to be interacting 
via the same type of epistasis almost exclusively. Most 
recently, there have been efforts in introducing measures 
to quantify monochromaticity [H, [23| . 
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In this work, we are motivated by these findings to ask 
how the monochromaticity feature emerged in nature. 
We want to gain insight into whether it can be explained 
as the outcome of a neutral evolutionary theory, driven by 
stochastic processes. To this end we devise network mod- 
els inspired by the duplication-divergence model [23427} 
that is shown to successfully capture the heterogeneous 
network properties and other important topological fea- 
tures of protein interaction networks. Focusing on the 
genomic level, we propose neutral evolutionary models in 
which epistatic networks of genes are constructed via sim- 
ilar duplication-divergence processes for the epistatic in- 
teraction signs. Signs of links are copied as genes are du- 
plicated and they are switched and/or rewired with finite 
probability. We generate ensembles of epistatic networks 
using these models and assess their monochromaticity 
using the PRISM algorithm [l2| which makes use of ag- 
glomerative clustering to detect modules. The number of 
monochromaticity violations Qmoduie of these networks 
are compared to their randomized versions where the 
topology is preserved and the epistatic signs are shuf- 
fled. We show that the Qmoduie distribution for epistatic 
network ensembles generated by these models are cen- 
tered around zero indicating monochromaticity, and are 
statistically signiflcantly different from their randomized 
counterparts, which are not monochromatic. In particu- 
lar, the last model that we propose in which we assume an 
epistasis schematically mediated by the interaction with 
the metabolic network, displays a clear separation of its 
conflict number distribution from that of the randomized 
ones and is strictly monochromatic. The degree distribu- 
tions of each of these models display a fat tailed behavior, 
with the degree distributions of positive and negative in- 
teractions showing no deviation from each other as there 
is no bias towards any sign in any of the models. 

This paper is structured as follows: In Section II we 
briefly revise the literature on duplication-divergence (D- 
D) models, summarize the PRISM algorithm and de- 
scribe in detail the D-D inspired evolutionary models that 
we use to construct the epistatic networks. We study the 
monochromaticity conflict distributions of ensembles of 
networks produced by these models and their randomized 
counterparts. We flnally make our concluding remarks in 
Section III. 



II. MONOCHROMATICITY IN 
DUPLICATION-DIVERGENCE INSPIRED 
EPISTATIC MODELS 

Gene dupHcation and the subsequent mutations (di- 
vergence) is believed to be one of the most crucial 
mechanisms driving evolution [28|. In the past decade, 
mathematical models that make use of the duplication- 
divergence mechanism have been proposed and studied 
analytically (23-[27j. These models were successful in 
capturing the basic topological properties of protein in- 
traction networks such as scale-free degree distribution. 



small- world properties, modular architecture and robust- 
ness against random node removal. The duplication parts 
of these models are carried out in a similar way to each 
other. First, a random gene is selected to be duphcated, 
which is called a target gene. The duplicated new gene 
then acquires all of the interactions of the target gene. 
This makes sense in the proteomic viewpoint as the in- 
teractions of proteins are determined by their structure, 
which remains unchanged after duplication. The diver- 
gence part, on the other hand, is model speciflc and 
shows some variation from model to model in terms of 
the exact divergence mechanism. Usually, what is under- 
stood by divergence is the removal of interactions with 
some flnite probability. Duplication-divergence models 
are divided into two classes according to the link re- 
moval mechanism. In asymmetric models [29| only 
the duplicated genes can lose links whereas in symmet- 
ric duplication-divergence models [25|, both the original 
and duplicated genes can lose their interactions. In some 
of the duplication-divergence models, this removal can be 
accompanied by rewiring of the links or formation of new 
links with some other flnite probability. 

In this paper we implement three epistatic network 
models based on the duplication-divergence mechanism 
with additional focus on the interaction signs. The 
epistatic interactions are mediated by metabolic and 
protein interaction networks. Therefore a duplication- 
divergence model for the evolution of epistatic interac- 
tions implements a minimal hypothesis on their dynam- 
ical evolution. We start from a simple initial condition 
common to all of the three models. During the duplica- 
tion phase, the interaction type is preserved whereas in 
the case of rewiring, the epistatic signs of the new links 
are chosen as detailed in the following subsections for 
each model. 

The monochromaticity - the feature of clusterability 
into modules that interact via the same type of epistasis 
- of each model is studied using the PRISM algorithm 
which was originally developed to study the yeast 
metabolic network. Belonging to a class of hierarchical 
clustering methods known as agglomerative clustering, 
the algorithm starts out with each element in a cluster 
of its own. It assigns a conflict score Cx^y to every clus- 
ter pair (x,?/), which is total number of "mixed sign" 
or nonmonochromatic links and then at each step deter- 
mines the set of cluster pairs with the minimum conflict 
score 

Cm = min{Cx,y} 

x,y 

It then selects from this set the one pair of clusters with 
the highest proximity (determined by the density of links 
between them) and combines them. This process is re- 
peated until the whole network is covered. Finally, a 
total monochromaticity violation number Qmoduie which 
is the sum of the minimum conflict scores over all the 
steps 

Qmoduie — / Cm 



is assigned to the final clustering solution. A total con- 
flict number of zero means a fully monochromatic solu- 
tion. The result of this analysis on the yeast metabolic 
network reported in [l2| is that the unsupervised organi- 
zation of the real yeast metabolic network yielded a fully 
monochromatic solution. Unsupervised in this context 
means having no prior knowledge of the genes' functional 
annotations. Moreover, this unsupervised organization of 
the epistatic network resulted in modules that agree well 
with existing gene annotations. A supervised version of 
this analysis was also carried out in this study for pre- 
assigned functional annotation groups, which also yielded 
a statistically significant enrichment of monochromatic 
interactions. 

In this paper we take a similar approach in which we 
investigate the distributions of the total conflict score 
Qmoduie for the networks constructed by duplication- 
divergence models and their randomized counterparts. 
The color randomization scheme consists of the pairwise 
flipping of signs so that the topology of the network is 
preserved as well as the total number of positive and neg- 
ative interactions. Following the unsupervised approach, 
we don't make any assumptions as to the function of 
genes. Since the duplication-divergence models are them- 
selves of stochastic nature, we represent their total con- 
flict score as a distribution as well. 



A. Model A 

The first model that we consider for the neutral evo- 
lution of epistatic networks is a simple generalization of 
the model introduced by Sole et al. [2^. The Sole model 
is one of the first models to be developed as a growing 
network model to simulate proteome evolution. It has 
two free parameters 5 and a which are the probability 
of removal of the links of the duplicated gene and the 
probability of the duplicated gene to form links with any 
of the previously unlinked genes, respectively. It was 
shown to successfully capture structural properties of 
the yeast proteome such as the degree distribution (U, 
sparse nature and modularity [30|. Here we generalize 
this model in order to capture a simple neutral evolution 
of the epistatic network. We start with simple initial 
conditions: at the beginning the network is initialized 
with three nodes and two links so that £^12 = 1 and 
£^13 = — 1 where E is the interaction matrix composed 
of the elements —1 (negative epistasis), 1 (positive 
epistasis) and (no interaction). The epistatic imple- 
mentation of this model for a network size of n consists 
of the iteration of the following set of rules (Fig. [1]). 



(i) duplication: A target node Vi \s chosen at ran- 
dom and the replicated node Vr = ^n+i acquires all 
the links {e^^j} of the target node where {vj} is the 
set of neighbors of Vi. The signs are copied during 
duplication so that Erj = Eij. 
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FIG. 1: Duplication-Divergence scheme of the Model A. The 
target node is indicated by the small arrow and the duplicate 
node is in grey. Green links denote positive interactions and 
red links negative interactions. Interactions are duplicated 
with the same signs; removed (dashed lines) with probability 
6 and linked (dotted lines) to previous nodes with random 
signs with probability a. 

(ii) divergence 1 (deletion): Each of the links Cr^j of the 
replicated node is deleted with probability 5. 

(iii) divergence 2 (addition of new links): The rephcated 
node Vr is connected with a random sign to all of 
the previously unconnected nodes with probability 
a, which is typically a small number. 

The parameters used in this work, namely (5 = 0.5 and 
a = 0.0002 were chosen so that they are compatible with 
the estimated average connectivity from mean-field cal- 
culations [24] and experimental data. 

An instance of the epistatic network generated by 
Model A is shown in Fig. [2l The total degree distribution 
and sign degree distributions (Fig. [3]) of this model both 
display fat tailed distributions with an exponential cutoff 
although it is easy to disrupt the power law distribution 
of this model by increasing the rewiring probability, in 
which case the distribution becomes exponential similar 
to that of a random growing network. Furthermore, both 
positive and negative interaction degree distributions fol- 
low the same trend as there is no bias in the model to- 
wards any sign. 




FIG. 2: A sample network given by the model A for (5 = 0.5 
and a = 0.0002. 

The distribution of the total number of conflicts of 
the networks generated by Model A and their random 
counterparts is shown in Fig. |4j To get this distri- 
bution, a total number of 100 networks were generated 
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FIG. 3: The total (left) and sign (right) degree distributions 
of the model A for ^ = 0.5 and a = 0.0002. 



and subsequently randomized. The resulting histogram 
shows a clear separation of the means in the duplication- 
divergence ensemble (shown in black) and randomized 
ensemble (shown in white). The conflict score of model A 
is very close to zero (0.91 ±1.45) whereas the correspond- 
ing randomized distribution has much larger spread. To 
test for the null hypothesis that the means of the distribu- 
tions are the same, or in other words the null assumption 
that the expected Qmoduie is given by the randomized 
network, we ran a statistical significance test on the dis- 
tributions, which yielded a p- value of 8.43 x 10~^^, which 
means that the distributions are statistically significantly 
different, with model A giving a monochromatic solution 
whereas the randomized version isn't. 
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FIG. 4: The total monochromaticity violation number 
Qmoduie distributions of model A networks (black) and their 
randomized counterparts (white) for 100 realizations of each. 
The parameters used are ^ — 0.5 and a = 0.0002. 



B. Model B 

In Model A the epistatic interactions of different signs 
remain strongly localized in different regions of the net- 
work as it appears clearly in the graphical representation 
of the network (Fig. [2]). In this section we provide a first 
modification of the model which has the effect of main- 
taining the fat-tail degree distribution while at the same 
time mixing further the epistatic interactions of different 
signs. Model B has two free parameters and takes into 



account the probability to connect to one of the second 
neighbors of the target node in the case of removal of the 
duplicated node. The duplication (i) step is the same as 
in the previous model whereas the divergence is different. 
Here in the following we give the precise definition of the 
algorithm we have implemented for Model B. The iter- 
ative steps are summarized in Figure [5l We start with 
the same initial conditions as in model A, i.e. we start 
with three nodes and two links such that the epistatic 
interactions are given by £^12 = 1 and £^13 = — 1 where 
E indicates the sign of the epistatic interaction. 



FIG. 5: Duplication-Divergence scheme of model B. The tar- 
get node is indicated by the small arrow and the duplicate 
node is in grey. Green links denote positive interactions and 
red links negative interactions. Interactions are duplicated 
with the same signs; removed (dashed lines) with probability 
p and in the case of removal, rewired to one of the second 
neighbors of the target node with probability q according to 
the sign convention described below. 



(i) duplication: A target node '^^ G V is chosen at ran- 
dom and the replicated node Vr = ^n+i acquires all 
the links {eij} of the target node where {vj} is the 
set of neighbors of Vi. The signs are copied during 
duplication so that Erj = Eij . 

(ii) divergence 1 (deletion): Each of the links Cr^j of the 
replicated node is deleted with probability p. 

(iii) divergence 2 (rewiring): The replicated node Vr is 
connected to one of the second neighbors of the tar- 
get node with probability q in the case of removal 
of the duplicated link. The sign of the rewired 
link to the second neighbor is the product of the 
signs of the first and second neighbors such that 



E, 



rk 



Ejj X E. 



This model suggests that when a gene is duplicated, 
newly emerging epistatic interactions follow a transitive 
sign rule (divergence 2 mechanism). Although epistatic 
interactions are known to be non-linear and in gen- 
eral a transitivity of the sign might not be the rule, it 
is possible that for epistatic interactions mediated by 
protein-interaction networks this mechanism for diver- 
gence might effectively take place. A sample network 
generated by Model B is given in (Fig. [6j). The free 
parameters for this model were chosen as p = 0.4 and 
q = 0.01 for removal and rewiring, respectively. Here the 
network has a more mixed topology due to the rewiring 
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to the second neighbors, therefore the monochromaticany 
clusterable architecture is not evident from the figure of 
the network. Still, the scale free degree distribution is 
preserved (Fig. [7j) as the rewiring is made with only one 
of the second neighbors of the target node, per removal 
of the duplicate link. Positive and negative degree distri- 
butions are again indiscernible. 




FIG. 6: A sample network given by model B for p = 0.4 and 
q = 0.01. 
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FIG. 7: The total (left) and sign (right) degree distributions 
of model B for p — 0.4 and q — 0.01. 

Duplication-divergence model B and randomized net- 
works have the Q module distribution as shown in (Fig. [8j). 
An ensemble of 100 networks was used to collect statis- 
tics on monochromaticity for this model. In Figure [8] the 
total conflict distribution of model B is presented and it 
appears as slightly more spread than Model A although 
it is still centered close to zero (3.83 ± 4.36). The corre- 
sponding randomized distribution is also more separated 
(61.41 ±53.22). The p- value under the same null hypoth- 
esis as the Model A is 1.26 x 10~^^. This also complies 
with the results above, indicating monochromatic sepa- 
rability of model B. 

C. Model C 

This model is a generalization of the basic model which 
preserves the fat-tail degree distribution and implies a 
well-mixed distribution of epistatic interactions with dif- 
ferent signs. The model is a very stylized one and it 
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FIG. 8: The total monochromaticity violation number 
Qmoduie distributions of model B networks (black) and their 
randomized counterparts (white) for 100 realizations of each. 
The parameters used are p = 0.4 and q — 0.01 

assumes that the metabolic network is the network un- 
derlying epistatic interactions. Therefore the fitness of 
the model organism is given by the biomass production 
during metabolism. In this metabolic network we have 
a flux of mass from the substrates to the production of 
biomass mediated by several metabolic reactions regu- 
lated by metabolic genes. A mutation on a gene can 
be seen as the knockout of a reaction perturbing the pro- 
duction of the biomass. A simple schematic argument as- 
sumes that buffering epistatic interactions between genes 
emerge when the genes regulate metabolic reactions that 
are "in series" in the metabolic network (Fig. [9|). Sim- 
ilarly aggravating epistatic interactions are supposed to 
be present between genes that regulate reactions that are 
"in parallel" in the metabolic network (Fig. [9]). In model 
C we consider evolution of epistatic interactions related 
to the metabolism. We assume that a duplicated gene 
might either regulate a reaction which is "in series" or 
"in parallel" with the reaction regulated by the original 
template gene. 




FIG. 9: Schematic view of the metabolic network as a network 
that sustains a flux of mass from the substrate to the pro- 
duction of biomass. Metabolic genes regulating the chemical 
reactions of the network have a buffering epistatic interaction 
if they are regulating chemical reactions which are "in series" 
in the metabolic network and have an aggravating interaction 
if they are regulating reactions that are "in parallel" in the 
metabolic network. 

We start with the same initial conditions as in model 




FIG. 10: Duplication-Divergence scheme of model C. The du- 
plicate node is in grey. Series (double green) and parallel 
(double red) reactions in the metabolic network are estab- 
Ushed between target and duplicate nodes with equal prob- 
ability. Signs of the duplicated links are chosen according 
to the definition of the model. Removed interactions (with 
probability p) are shown by dashed lines. 

A and B, i.e. we start with three nodes and two links 
such that the epistatic interactions are given by £^12 = 1 
and £^13 = — 1 where E indicates the sign of the epistatic 
interaction. At each iteration we perform the subsequent 
steps (Fig. flQ)) : 

(i) duplication (with series or parallel connections in 
the metabolic network): A new duplicated node 
Vr = Vn-\-i is "in series" {Eri = 1) or "in parallel" 
(Eri = — 1) with respect to a randomly chosen tem- 
plate node Vi ^ V in the metabolic network with 
equal probability. If the connection is "in series", 
the replicated node acquires all the links {eij} of 
the target node with the same sign, i.e. E^j = Eij. 
If the connection is "in parallel" , the replicated node 
acquires all the links {eij} of the target node with 
the opposite sign, i.e. E^j = —Eij. 

(ii) divergence (deletion): Each of the links Crj of the 
replicated node is deleted with probability p. 

AA 

FIG. 11: The two triangular motifs generated by model C. 

A typical network constructed by Model C is shown in 
(Fig. [T2]) . The only free parameter is the probability 
of removal of the duplicated links, which is chosen to 
be p = 0.6. The epistatic signs are thoroughly mixed 
in this network, due to the equally probable assignment 
of positive and negative interactions between the target 
node and the duplicated node, and the ensuing sign con- 
vention of the duplicated links. This model permits the 
generation of two types of triangular motifs, namely pnn 
and ppp triangles (Fig. [TT]) where p and n indicate re- 
spectively positive (buffering) and negative (aggravating) 
interactions. The degree distributions of this model are 
fat tailed as well despite the presence of some skewness 
in the power law distribution. The degree distributions 
of interactions with different signs follow the same trend 

(Fig. m- 
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FIG. 12: A sample network given by model C for p = 0.6. 
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FIG. 13: The total (left) and sign (right) degree distributions 
of model C for p = 0.6. 



An ensemble of 150 networks was used for the total 
conflict distribution. Compared to the other two models, 
this model is much more defined in terms of monochro- 
maticity and gives a fully monochromatic solution (Fig. 
[T4|) . The separation of the duplication divergence model 
from its randomized counterpart (80.16 ± 18.80) is ap- 
parent. 
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FIG. 14: The total monochromaticity violation number 
Qmoduie distributions of model C networks (black) and their 
randomized counterparts (white) for 150 realizations of each. 
The parameter is chosen as p = 0.6 



III. CONCLUSION 

In conlusion in this paper we have investigated the 
modular structure and monochromatic clusterability of 
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epistatic networks generated using stochastic models in- 
spired by the duphcation-divergence method in order 
to test for a neutral evolutionary hypothesis. We have 
shown that these genetic interaction networks have to- 
tal monochromatic violation score distributions that are 
close to zero and are statistically significantly different 
from their randomized counterparts. This fact implies 
that these networks comply with the monochromaticity 
feature observed in real genetic interaction data of model 
organisms such as the budding yeast. 

Overall, the investigation of models based on 
the duplication-divergence mechanism is a good first 
step toward understanding the evolutionary origins of 
monochromaticity although it is important to stress that 
these mathematical models are very simplistic and fur- 
ther investigation into real biological networks is called 
for to extract information about what features are more 
relevant to monochromaticity in order to be able to tweak 
the details of these models. The parameter space of the 
three proposed models were explored before the Qmoduie 
distributions were plotted and the free parameters were 
chosen so as to reproduce the general topological proper- 
ties (scale-free degree distribution, average connectivity) 
of model organisms. Therefore it is fair to say that the 
results presented in this study constitute a snapshot of 
the models since modularity comes about in a relatively 
narrow range of parameters. Duplication-divergence does 
not account for the preference of these specific parame- 



ters. It is likely that they might be the result of some 
selection for sparse graphs in nature [3l|. So biological 
modularity cannot be explained solely by natural selec- 
tion or neutral evolution but rather as the result of con- 
tributions from both - an example being neutral mech- 
anisms such as duplication-divergence whose parameters 
are tweaked by natural selection. The same is most likely 
valid for monochromaticity. In fact it is difficult to say 
that neutral evolution is the only explanation underly- 
ing it, although our findings with duplication-divergence 
models do corroborate its role. 

One possible future direction to take, as far as assess- 
ing what the relevant structural characteristics in real bi- 
ological networks are, would be to do a motif significance 
investigation on real data. Motifs are highly represented 
subgraphs in networks that occur statistically more fre- 
quently than they should in the randomized counterparts 
of networks [32] . Although the immediate functional sig- 
nificance of individual motifs has been a disputed topic, 
it is still informative about the network structure. In- 
teresting parallels might be drawn between the motifs 
of real networks and implemented mathematical models. 
Model C in particular calls for such an analysis for dif- 
ferent triangle types (4 overall). A comparison between 
these particular motifs in the yeast interaction network 
and networks generated by Model C can help improve 
the already promising ability of this model to reproduce 
monochromaticity in a biologically motivated way. 
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